clear;
clc;

data = [
36.21396127,120.08041820;
36.21384291,120.08062440;
36.21383899,120.08092840;
36.21384037,120.08121860;
36.21384005,120.08154620;
36.21384199,120.08180660;
36.21384278,120.08213660;
36.21392075,120.08225150;
36.21409733,120.08226650;
36.21434936,120.08225520;
36.21443832,120.08215390;
36.21444200,120.08192100;
36.21444049,120.08162530;
36.21444216,120.08132740;
36.21443959,120.08097680;
36.21443931,120.08067080;
36.21441290,120.08048700;
36.21430002,120.08041290;
36.21409796,120.08041360;
]; 

wgs_lon = data(:,2);
wgs_lat = data(:,1);

[gcj_lon, gcj_lat] = wgs84_to_gcj02_batch(wgs_lon, wgs_lat);

data_gcj = [gcj_lat,gcj_lon];

fprintf('%.8f,%.8f,',data_gcj');

figure;
subplot(1, 1, 1);
plot(gcj_lat(:), gcj_lon(:), 'mo--', 'DisplayName', '30个采样点'); % 均匀采样点
title('程序验证生成点');
xlabel('X');
ylabel('Y');


function [gcj_lons, gcj_lats] = wgs84_to_gcj02_batch(wgs_lons, wgs_lats)
% WGS84坐标批量转GCJ02坐标（火星坐标系）
% 输入: 
%   wgs_lons - WGS84经度数组（列向量或行向量）
%   wgs_lats - WGS84纬度数组（列向量或行向量）
% 输出: 
%   gcj_lons - GCJ02经度数组
%   gcj_lats - GCJ02纬度数组
%
% 示例:
%   wgs_lons = [120.0804; 116.404; 121.4737];
%   wgs_lats = [36.2142; 39.915; 31.2304];
%   [gcj_lons, gcj_lats] = wgs84_to_gcj02_batch(wgs_lons, wgs_lats);

    % 检查输入数据是否长度一致
    if length(wgs_lons) ~= length(wgs_lats)
        error('经度和纬度数组长度不一致！');
    end

    % 初始化输出数组
    num_points = length(wgs_lons);
    gcj_lons = zeros(num_points, 1);
    gcj_lats = zeros(num_points, 1);

    % 逐点转换
    for i = 1:num_points
        [gcj_lons(i), gcj_lats(i)] = wgs84_to_gcj02_single(wgs_lons(i), wgs_lats(i));
    end
end

function [gcj_lon, gcj_lat] = wgs84_to_gcj02_single(wgs_lon, wgs_lat)
% 单点WGS84转GCJ02
    % Krasovsky 1940椭球参数
    a = 6378245.0; % 长半轴
    ee = 0.00669342162296594323; % 扁率

    % 检查是否在中国境内
    if ~is_in_china(wgs_lon, wgs_lat)
        gcj_lon = wgs_lon;
        gcj_lat = wgs_lat;
        return;
    end

    % 计算偏移
    dlat = transform_lat(wgs_lon - 105.0, wgs_lat - 35.0);
    dlon = transform_lon(wgs_lon - 105.0, wgs_lat - 35.0);

    rad_lat = wgs_lat / 180.0 * pi;
    magic = sin(rad_lat);
    magic = 1 - ee * magic * magic;
    sqrt_magic = sqrt(magic);

    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrt_magic) * pi);
    dlon = (dlon * 180.0) / (a / sqrt_magic * cos(rad_lat) * pi);

    gcj_lat = wgs_lat + dlat;
    gcj_lon = wgs_lon + dlon;
end

function delta = transform_lat(x, y)
    % 纬度偏移计算
    ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x));
    ret = ret + (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret = ret + (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0;
    ret = ret + (160.0 * sin(y / 12.0 * pi) + 320.0 * sin(y * pi / 30.0)) * 2.0 / 3.0;
    delta = ret;
end

function delta = transform_lon(x, y)
    % 经度偏移计算
    ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x));
    ret = ret + (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret = ret + (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0;
    ret = ret + (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0;
    delta = ret;
end

function in_china = is_in_china(lon, lat)
    % 判断坐标是否在中国境内（粗略边界）
    in_china = (lon >= 72.004 && lon <= 137.8347 && lat >= 0.8293 && lat <= 55.8271);
end